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Abstract 

The long-known, widely documented inverse relationship between body size 
and environmental temperature (“temperature-size rule”) has recently led to 
predictions of body size decline following current climatic warming (“size 
shrinking effect”). For keystone pollinators such as wild bees, body shrinking 
in response to warming can have significant effects on pollination processes 
but there is still little direct evidence of the phenomenon because adequate 
tests require controlling for confounding factors linked to climate change 
(e.g., habitat change). This paper assesses the shrinking effect in a community 
of solitary bees from well-preserved habitats in the core of a large nature 
reserve experiencing climatic warming without disturbances or habitat 
changes. Long-term variation in mean body mass was evaluated using data 
from 1704 individual bees (137 species, 27 genera, 6 families) sampled over 
1990-2023. Climate warmed at a fast rate during this period, annual mean of 
daily maximum temperature increasing 0.069°C/year on average during 
2000-2020. Changes in bee body mass verified expectations from the size 
shrinking effect. The mean individual body mass of the community of solitary 
bees declined significantly, irrespective of whether the analysis referred to the 
full species sample or only to the subset of species that were sampled in both 
the old (1990-1997) and recent (2022-2023) periods. On average, body mass 
declined ~0.7%-year™', or an estimated average cumulative reduction of 
~20 mg per individual bee from 1990 to 2023. Proportional size reduction was 
greatest among large-bodied species, ranging from around —0.6%-year~’ for 
the smallest species to —0.9%-year~* for the largest ones. Declining rate was 
steeper for cavity-nesting than ground-nesting species. The pollination and 
mating systems of bee-pollinated plants in the study region are probably 
undergoing significant alterations as a consequence of supra-annual decline in 
bee body mass. 
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INTRODUCTION 


Intra- and interspecific variation in individual body size 
is often predictably related to environmental tempera- 
ture. Early recognition of such relationships led to the 
formulation of some classical ecological rules such as 
Bergmann’s, Allen’s, or the “temperature-size rule” 
(Atkinson, 1994; Ray, 1960). The temperature-size rule, a 
“widely documented and poorly understood pattern” 
(Klok & Harrison, 2013), predicts an inverse relationship 
across conspecific individuals between body size and the 
environmental temperature experienced during develop- 
ment (Atkinson, 1994; Forster & Hirst, 2012; Klok & 
Harrison, 2013). While the proximate and ultimate causes 
of the temperature-size rule are still far from being fully 
elucidated (Verberk et al., 2021; Walters & Hassall, 2006), 
it is a matter of fact that the rule has proven true in the 
majority of ectothermic organisms where it has been 
investigated (Atkinson, 1994; Klok & Harrison, 2013). 
Such a predictable relationship between body size and 
environmental temperature has recently acquired partic- 
ular relevance in current ecological scenarios affected by 
climate warming, motivating predictions of reduced body 
size in response to rising temperatures (“size shrinking 
effect”; Ohlberger, 2013; Sheridan & Bickford, 2011; 
Verberk et al., 2021). 

Given the many ecological and evolutionary conse- 
quences of body size (Bonner, 2006), its predicted decline 
due to climatic warming can trigger a complex cascade of 
effects whose details will depend on the ecological func- 
tionality of the organisms involved. In the case of key- 
stone pollinators such as wild bees, body shrinking can 
alter crucial aspects of the pollination process such as for- 
aging range, pollen load size and diversity, and pollen 
carryover, all of which depend on body size (Cullen et al., 
2021; Földesi et al., 2021; Greenleaf et al., 2007). Limited 
evidence suggests that recent changes in bee sizes con- 
form to expectations from the temperature-size rule, but 
studies conducted so far refer to few species, consider lin- 
ear measurements as a proxy for body size rather than 
actual body mass, and/or the putative effects of tempera- 
ture can be confounded with those of other factors corre- 
lated with climate change such as urbanization, land use 
changes, or habitat fragmentation (Barrett & Johnson, 
2023; Garlin et al., 2022; Grab et al., 2019; Kelemen & 
Rehan, 2021; Oliveira et al., 2016; Suni & Dela Cruz, 
2021). In this paper, we assess the size shrinking effect 


expected from the temperature-size rule by examining 
long-term changes in individual body mass in a diverse 
community of solitary bees sampled over three decades 
in Mediterranean habitats from a large protected area. 
Our study region is undergoing climatic warming but not 
disturbances or habitat changes that could confound the 
effects of the warming climate. 


MATERIALS AND METHODS 


As part of other studies (Herrera, 1995, 1997; Herrera et al., 
2023; C. M. Herrera, unpublished data), a total of 1704 wild 
bees from 137 species, 27 genera, and six families 
(Andrenidae, Apidae, Colletidae, Halictidae, Megachilidae, 
Melittidae; see Appendix S1: Table S1, for species and sam- 
ple sizes) were hand-netted in the field during 1990-1997 
(“old period” hereafter; N = 473 individuals, 47 species) 
and 2022-2023 (“recent period”; N = 1231 individuals, 
130 species), in 55 localities from the core area of the 
Sierras de Cazorla, Segura y Las Villas Natural Park, a 
2090 km? protected area in Jaén Province, southeastern 
Spain. A map of sampling locations is shown in 
Appendix S2: Figure S1. The sampled area is characterized 
by well-preserved habitats and outstanding biological diver- 
sity (GOmez Mercado, 2011; Médail & Diadema, 2009; 
Molina-Venegas et al., 2015). Our bee sampling sites were 
located at elevations between 740 and 1988 m, and did not 
undergo any obvious disturbance or habitat change over 
the 30-year period considered here which could have 
affected wild bee populations (e.g., wildfires, introduction 
of managed honeybee colonies, arrival or expansion of 
invasive plants). Total annual precipitation fluctuated 
widely over the study period but without a significant tem- 
poral trend (Herrera, 2019). 

Netted bees were placed individually into sealed 
microcentrifuge tubes, kept in the dark in a portable 
refrigerator at 4-8°C, and brought to the laboratory 
within a few hours, where they were weighed to the 
nearest 0.1 mg always using the same Sartorius 1602MP8 
analytical balance. The species sampled represented 
about one third of all species of bees occurring in the 
region (Ortiz-Sanchez et al., In press) and encompassed 
the whole range of body sizes (range of species 
means = 6-544 mg). Most bees weighed belonged to the 
genera Andrena (864 individuals, 40 species), Anthophora 
(167, 13), Colletes (139, 4), Osmia (112, 6), Xylocopa 
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(90, 3), and Anthidium (64, 2) (Appendix S1: Table S1). 
All bee species sampled were solitary, with the minor 
exception of a small proportion of data from primitively 
eusocial Halictids in the genera Lasioglossum and 
Halictus (15 individuals from five species; sociality data 
for Halictidae based on Appendix A of Gibbs et al., 2012). 
Each species was assigned to either ground-nesting (“soil 
excavators” sensu Danforth et al., 2019) or cavity-nesting 
(comprising “wood excavators,” “renters” and “above- 
ground builders” types of Danforth et al., 2019) catego- 
ries, using information from Michener (2000) and 
Danforth et al. (2019). Brood parasites (Coelioxys, 
Epeolus, Melecta, Nomada, Thyreus) were assigned to the 
same nesting class as their most frequent hosts. Body 
mass data of all individual bees included in this study, 
along with associated metadata (sex, nesting category, 
sampling date and site, and elevation and geographical 
coordinates of sampling site), are publicly available 
(Herrera, 2023). 

Recent increase in environmental temperature has been 
well documented for the southeastern Iberian Peninsula 
(Fernandez-Montes & Rodrigo, 2015; Gonzalez-Hidalgo 
et al., 2015). To corroborate this trend at the spatial and 
temporal scale of this study, daily maximum and minimum 
temperature data for 2000-2022 were gathered for 
10 weather stations 15-60 km away from our bee sampling 
sites (Red de Información Agroclimatica de Andalucia; 
https://www.juntadeandalucia.es/agriculturaypesca/ifapa/ 
riaweb/web/; last accessed 23 October 2022). Relevant 
details and location map of weather stations are shown 
in Appendix S2: Figure S1, Table S1. 


Data analysis 


A similar analytical framework based on fitting linear 
mixed effects models to the data (“mixed models” hereaf- 
ter; Bolker, 2015; Harrison et al., 2018) was adopted to 
test for supra-annual trends in both environmental tem- 
perature and body mass of solitary bees. In the case of 
temperatures, random intercepts models were fitted to 
daily maximum and daily minimum temperature data for 
every weather station. These models had year as the sin- 
gle fixed effect (treated as a continuous, numerical vari- 
able) and Julian date (=days since 1 January) as random 
effect. To circumvent the annual periodicity underlying 
weather data, Julian date was treated as an unordered 
random factor with 365 levels. Fixed-effect parameters 
obtained from these models provided estimates of the 
rate of change of annual means for daily maximum and 
daily minimum temperatures. For bee body mass, several 
random slopes-random intercepts models were fitted to 
individual mass data (logio transformed), as detailed in 


the next paragraph. All these models had year of capture 
(centered and scaled to facilitate model convergence and 
interpretation of fixed effects; Harrison et al., 2018) as a 
fixed effect predictor, and bee species as random effect. 
Since sexual size dimorphism is frequent in solitary bees 
(Danforth et al., 2019), sex and its interaction with year 
of capture were included in each model as fixed-effect 
covariates. Preliminary analyses including Julian date 
and altitude of sampling site as additional fixed-effect 
covariates did not improve the fit of any model signifi- 
cantly. For simplicity, these two variables were omitted 
from the analyses reported here. 

Most bee species were sampled on only one or a few 
years, which resulted in a sparse species x year data 
matrix (Appendix S1: Table S1). This was mainly a conse- 
quence of enhanced sampling effort in the recent period, 
but it could also conceal possible long-term change in 
bee community composition in response to climatic 
warming (e.g., increased representation of small-sized 
species and/or rarefaction of large-sized ones; Herrera, 
2019; Herrera et al., 2023). Although linear mixed models 
are inherently robust to data sparseness (i.e., sparsely 
sampled levels of random effects, Bolker, 2015) and viola- 
tion of distributional assumptions (Schielzeth et al., 
2020), long-term changes in bee community composition 
could still bias our analyses because of changes in the 
sampling universe, and hence inference space, between 
the old and recent periods. To evaluate this possibility, 
mixed models were fitted to data sets that differed in tem- 
poral scope and bee species composition. The two main 
analyses considered the complete temporal scope of our 
study (1990-2023), but differed in the species composi- 
tion of the sample. In one case, the model was fitted to 
all the data (N = 1704, 137 species), while in the other, 
only species that were sampled in both the old and recent 
periods were included (N = 1385, 40 species). Two addi- 
tional, supplementary analyses were conducted which 
used data for all bee species but considered distinct tem- 
poral scopes: one was fitted to data from the old period 
(1990-1997; N = 473, 47 species) and the other to data 
from the recent one (2022-2023; N = 1231, 130 species). 

All statistical analyses were carried out using the R 
environment (R Core Team, 2022). Mixed models were 
fitted with the Imer function in the Ime4 package (Bates 
et al., 2015). Model fit adequacy was assessed using the 
check_model function in the performance package 
(Liidecke et al., 2021). In the analyses of weather data, con- 
fidence intervals of fixed-effect parameters were obtained 
by bootstrapping (function bootstrap_model from the 
parameters package; Ltidecke et al., 2020), and p values for 
tests of the same hypothesis for multiple weather stations 
were corrected using the Benjamini-Hochberg procedure 
(p.adjust function, stats package). Statistical significance of 
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fixed effects was assessed using Wald y? tests (Anova 
function, car package; Fox & Weisberg, 2019). The function 
ggpredict from the ggeffects package (Lüdecke, 2018) was 
used to compute marginal effects of year on mean tempera- 
ture and mean bee body mass. Random slopes from the 
mixed model fitted to all bee mass data, which reflected 
the variation among levels of the random effect (bee spe- 
cies) in the effect of the predictor (year) on the response 
variable (body mass), were obtained using the coef func- 
tion. These species-specific declining rates were then 
related to interspecific variation in body mass and nesting 
habit. 


RESULTS 


The study area underwent significant warming over the 
period 2000-2022. This trend was mainly due to a steep 
increase in yearly means of daily maximum temperatures, 
which took place consistently at all weather stations 
(Figure 1). All mixed models fitted to daily maximum 
temperature data yielded statistically significant, positive 
temperature/year relationships, with model-estimated 
slopes averaging +0.069°C/year and ranging between 
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FIGURE 1 Linear trends in annual means for daily maximum 
and daily minimum temperatures for 10 weather stations near and 
around the bee sampling sites (see Appendix S2: Figure S1 for a 
map). Each line is the prediction obtained from the linear mixed 
model fitted to daily temperature data for one station, with year as 
the single fixed effect and Julian date (unordered, qualitative factor) 
as random effect. A summary of analytical results is given in 
Appendix S2: Table S1. 


+0.042 and +0.12°C/year depending on the station 
(see Appendix S2: Table S1, for parameter estimates and 
confidence intervals). Increases in mean daily minimum 
temperature also took place, but were less marked 
(Figure 1). Statistically significant relationships between 
daily minimum temperature and year occurred in only 
six weather stations, and in these cases, the estimated 
slopes averaged +0.036°C/year and ranged between 
+0.012 and +0.088°C/year (Appendix S2: Table S1). 

Mean individual body mass of the community of soli- 
tary bees declined significantly over the 1990-2023 sam- 
pling period, irrespective of whether the analysis referred 
to all species or only to those that were sampled in both 
the old and recent period (Table 1). The two main ana- 
lyses revealed a negative effect of year on body mass after 
Statistically accounting for significant sexual size dimor- 
phism and allowing for interspecific variation in inter- 
cepts and slopes of the body mass/year relationship 
(Table 1). The two models yielded remarkably similar 
parameter estimates, and fitted well the (logio 
transformed) individual body mass data and model 
assumptions (Appendix S3: Figures S1 and S2). The two 
supplementary analyses, fitting separate mixed models 
on data from the old and recent sampling periods, also 
revealed statistically significant trends of declining body 
mass within each period (Table 1). Results for the 
sex X year interaction effect were inconsistent among 
models, being statistically nonsignificant when the data 
encompassed the whole study period and significant in 
each of the two within-period analyses (Table 1). Only 
the main effect of year on body mass estimated from the 
two models fitted to data from the 1990 to 2023 period 
will be considered hereafter. 

Because of the transformations applied to year and 
mass, the estimated slopes of the body mass/year rela- 
tionship over the 1990-2023 period for the whole data set 
and for the subset of species sampled in both periods 
(—0.0386 and —0.0378, respectively, Table 1) denote 
changes in log;o(mass) per year standard deviation unit 
and reflect multiplicative rather than additive changes. 
Back-transforming these slopes to the original scales of 
measurement yielded average body mass declines 
of 0.681%-year~' and 0.666%-year', for all the data and 
for species sampled in the two periods, respectively, or an 
estimated average cumulative reduction of ~20 mg per 
individual bee from 1990 to 2023 (Figure 2). 
Within-species decline in mean body mass for every 
species X sex combinations which was sampled in both 
the old and recent periods are depicted in Figure 3. 

Bee species differed in annual rates of body mass 
decline, as shown graphically in Figure 3 and analytically 
by the significant difference between mixed models with 
and without random slopes, both for the models fitted to 


ECOLOGY 5 of 9 


TABLE 1 Results of linear mixed models fitted to body mass data for individual bees sampled over 1990-2023. 


Analysis and data set Fixed effect Parameter estimate [SE] x p-value 
Main analyses (all study period, 1990-2023) 
All data (N = 1704, 137 species) Year —0.0385 [0.0056] 47.2 6.3e-12 
Sex (male) —0.3015 [0.0059] 2614.3 <2.2e-16 
Year x sex (male) 0.0116 [0.0060] 3.7 0.055 
Species sampled in old and recent periods Year —0.0377 [0.0056] 48.6 3.2e-12 
(N = 1385, 40 species) Sex (male) —0.3091 [0.0062] 2480.7  <2.2e-16 
Year x sex (male) 0.0077 [0.0061] 1.6 0.21 
Supplementary analyses (within periods) 
Old period, 1990-1997, all species Year —0.0208 [0.1627] 5.3 0.021 
(N = 473, 47 species) Sex (male) —1.4564 [0.2281] 795.3  <2.2e-16 
Year x sex (male) —0.7305 [0.1472] 24.6 6.9e-07 
Recent period, 2022-2023, all species Year —0.0302 [0.1177] 3.9 0.049 
(N = 1231, 130 species) Sex (male) 0.0376 [0.1122] 1658.9 <2.2e-16 
Year x sex (male) —0.5094 [0.1746] 8.5 0.0035 


Note: Model parameters are expressed in the transformed scales (body mass log), transformed, years scaled to mean = 0 and standard deviation = 1). 
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FIGURE 2 Mean estimated marginal effect of year on bee 
body mass, as predicted from random slopes-random intercepts 
mixed models with body mass as response variable, year, sex, and 
their interaction as fixed effects, and bee species as random effect 
(Table 1). Separate analyses were done on the whole sample (green) 
and on the subset of species which were sampled on both the old 
(1990-1997) and recent (2022-2023) periods (red). The original 
analyses were conducted on transformed variables, and the mean 
marginal effects shown here are the back-transforms to the original 
measurement scales. 


all data (x? = 23.3, df = 2, p = 8.8e-06; likelihood ratio 
test) and only to data from species sampled in the two 
periods (x° = 23.3, df=2, p= 8.8e-06; likelihood 


ratio test). For all data, species-specific estimated slopes 
ranged between —0.078 and + 0.007 log; (mass)/year 
standard deviation, and were negative in 136 out of the 
137 instances (99.3%). Body mass was a significant pre- 
dictor of interspecific variation in declining slope for the 
whole set of species (F135 = 28.9, p = 3.2e-07), but not 
for the subset of species sampled in both periods 
(F\.33 = 3.6, p = 0.065). In the first case, the heavier a bee 
species, the steeper the declining rate of body mass over 
the study period (Figure 4). Back-transforming the 
values plotted in Figure 4 to the original measurement 
scales, species-specific proportional declining rates 
ranged from around —0.6%-year~' for the smallest species 
to —0.9%-year~' for the largest ones. The relationship 
between species-specific declining slope and nesting 
habit (ground vs. cavity-nesting) was not statistically sig- 
nificant in the whole data set (F135 = 2.7, p = 0.10), but 
it was in the subset of species sampled in both periods 
(Fi 33 = 8.5, p = 0.006). In this latter data set, average 
declining rate was steeper for cavity-nesting than 
ground-nesting species (mean + SE = —0.0524 + 0.0046 
and —0.0374 + 0.0023, respectively). 


DISCUSSION 


The climate of our study area warmed significantly 
during the past few decades, as shown by increasing 
mean daily maximum and, to a lesser extent, minimum 
temperatures, thus confirming the general trend for the 
southeastern Iberian Peninsula (see references in the 
Introduction section) at the reduced spatial scale of our 
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FIGURE 3 Within-species paired comparisons of mean body mass for all sex x species combinations (N = 51, involving 37 species in 


11 genera) that were sampled in both the old (1990-1997) and recent (2022-2023) periods. Each line connects the old and recent averages for 


the same species. 


0.00 
oO 

2 -0.02 
O 
O 
2 

œ —0.04 
= 
O 
LL 

—0.06 

—0.08 

10 30 100 300 
Mean body mass (mg) 
FIGURE 4 Inverse relationship between the slopes of 


logio(body mass)/year relationships estimated from the mixed 
model fitted to all the data (Table 1) and the mean body mass of 
individual bee species (N = 137). Line is the least-squares 
regression. 


study. In agreement with expectations from the size 
shrinking effect, the warming climatic trend was con- 
comitant with a decline in the mean individual body 
mass of the regional community of solitary bees over 


the whole study period (1990-2023), irrespective of 
whether the whole data set or only the subset of spe- 
cies sampled in both periods were considered. The 
declining trend was even discernible within each of the 
old and recent periods, despite the decreased statistical 
power to detect trends due to reduction in sample sizes 
and, particularly, temporal ranges, which in the recent 
period consisted of just 2 years. Most remarkably, the 
four log;,(mass)/year fixed-effect estimates fell within a 
narrow interval, irrespective of temporal span and 
species composition of the data set (from —0.021 to 
—0.038, Table 1). Taken together, these findings allow 
to confidently rule out the possibility that the observed 
trend in body size reduction is a spurious effect of het- 
erogeneous sampling and possible long-term changes 
in bee community composition. We thus conclude that 
long-term reduction in mean community-level bee 
body mass was chiefly or exclusively a consequence of 
the pervasive reduction in body size experienced by 
individual bee species. 

Bees for this study were sampled in a well-preserved, 
protected area located >10 km away from urbanized or 
agricultural land; hence, observed body size reduction 
can be parsimoniously attributed to the effects of climate 
warming. These results also provide unique evidence to 
date of size shrinking in a diverse wild bee community 
based on body mass data rather than size metrics based 
on linear measurements of body parts, which are poor 
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predictors of intraspecific variation in body mass (Kendall 
et al., 2019) and can produce biased results because of 
heterogeneous responses of body parts to temperature 
(Klok & Harrison, 2013; Ray, 1960). Adding robustness to 
our results is the treatment of the different bee spe- 
cies as levels of a random effect. Mixed models allow 
to make inferences that apply to different populations 
of effects, or “inference spaces” (Schabenberger & 
Pierce, 2002). In the context of this study, the whole 
regional community of solitary bees represents the 
“broad inference space” and our conclusions refer 
specifically to that space, not just the set of species 
sampled. This means that model parameter estimates 
for fixed effects refer to temporal changes in mean 
individual body mass for the solitary bee community 
as a whole, not just the 137 species sampled (see 
Herrera, 2019, for further discussion on the value of 
treating species as random effect levels when investi- 
gating community-level trends). 

As it often happens with natural patterns conforming 
to the temperature-size rule (Klok & Harrison, 2013; 
Verberk et al., 2021), the possible mechanism(s) responsi- 
ble for the decline in body mass of solitary bees 
documented in this paper can only be tentatively 
suggested. The fact that estimated slopes for individual 
species were consistently negative suggests that the ulti- 
mate cause of size reduction was universal enough as to 
affect similarly to all species. Experimental studies under 
controlled conditions have documented inverse effects on 
body size of the temperature experienced during larval 
growth for some genera or species of bees included in this 
study (e.g., Osmia bicornis, Lasioglossum; Kamm, 1974; 
Kierat et al., 2017; Radmacher & Strohm, 2010). These 
findings are in line with the size-temperature rule, and 
point to the ubiquitous direct effects of rising tempera- 
tures as the principal cause of body shrinking of solitary 
bees in our study region. This interpretation is supported 
by our finding that, for species sampled in both periods, 
species-specific declining slopes were steeper on average 
for species nesting in above-ground cavities than for 
those nesting underground, since developing larvae in 
underground nests are expected to be better buffered 
against rising ambient temperatures than those in 
above-ground cavities owing to the insulatory properties 
of soil (Antoine & Forrest, 2021). In addition to its direct 
effect on larval development, however, long-term 
increase in ambient temperature may also have indirect 
effects on adult body size through effects on food supply. 
Climate warming can reduce the floral food resources 
available to bees (Moss & Evans, 2022; Takkis et al., 
2015), and impairment in quantity and quality of larval 
food provisions will also have a negative effect on adult 
body sizes (Chole et al., 2019). In Mediterranean-climate 


areas, supra-annual variation in rainfall can influence 
primary production and availability of floral resources, 
which could eventually translate into variations in bee 
body size. The influence of this effect on the patterns 
revealed by this study, if any, is probably negligible, 
because (1) no significant supra-annual trend in total 
annual rainfall is currently taking place in the study 
region, and (2) there is a significant trend towards pro- 
portionally more rainfall falling in January-June 
(Herrera, 2019; C. M. Herrera, unpublished data), which 
includes the period of year when precipitation is most 
influential on primary production in Mediterranean eco- 
systems (Bartsch et al., 2020). 

Larger bee species experienced the largest propor- 
tional reductions in body mass (see Oliveira et al., 2016, 
for similar patterns). Regardless of the factors accounting 
for this differential response, which cannot be addressed 
with the data available, size-dependent reduction in bee 
body size has two important implications. First, the close 
relationship linking body size and fecundity in insects 
(Honek, 1993) suggests that per-capita fecundity of the 
largest species of solitary bees may have declined signifi- 
cantly in our study region over the past few decades. This 
effect could account for the recent rarefaction of some 
large-sized species (e.g., Andrena assimilis, A. thoracica, 
Xylocopa cantabrita; C. M. Herrera, unpublished data). 
And second, larger species tend to forage over wider 
areas, perform more and more effective pollinations, and 
deposit larger and more diverse pollen loads with greater 
carryover (Cullen et al., 2021; Földesi et al., 2021; 
Greenleaf et al., 2007; Herrera, 1987). The fact that larger 
bees experience the steepest declines in body mass thus 
suggests that the pollination and mating systems of many 
bee-pollinated plants of our study region may currently 
be undergoing significant changes and, possibly, also 
some shift in selective regime. These effects are predicted 
to occur most frequently among plants with large flowers 
with restrictive corollas that are predominantly or exclu- 
sively pollinated by large bees (e.g., Antirhinum, Digitalis, 
Helleborus, Linaria; Herrera, 2020). 
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